Here are the details of the Bayesian multi-level modelling. Our general approach is:
Define Priors
In this section we will specify some priors. We then then use a prior-predictive check to assess whether the prior is reasonable or not (i.e., on the same order of magnitude as our measurements).
Fixed Effects
Our independent variable is a categorical factor with 16 levels. We will drop the intercept from our model and instead fit a coefficent for each factor level (\(y \sim x - 0\)). As our dependant variable will be log-transformed, we can use the priors below:
prior <- c(
set_prior("normal(0,2)", class = "b"),
set_prior("cauchy(0,2)", class = "sigma"))
Group-level Effects
We will keep the default weakly informative priors for the group-level (‘random’) effects. From the brms documentation:
[…] restricted to be non-negative and, by default, have a half student-t prior with 3 degrees of freedom and a scale parameter that depends on the standard deviation of the response after applying the link function. Minimally, the scale parameter is 10. This prior is used (a) to be only very weakly informative in order to influence results as few as possible, while (b) providing at least some regularization to considerably improve convergence and sampling efficiency.
Prior Predictive Check
Now we can specify our Bayesian multi-level model and priors. Note that as we are using sample_prior = 'only', the model will not learn anything from our data.
m_prior <- brm(data = d_eeg,
rms ~ wg-1 + (1|subject),
family = "lognormal",
prior = prior,
iter = n_iter ,
sample_prior = 'only')
We can use this model to generate data.
## Observations: 320,000
## Variables: 2
## $ key <chr> "P2", "P2", "P2", "P2", "P2", "P2", "P2", "P2", "P2", "P2"…
## $ value <dbl> 0.58997505, 2.38399527, 0.16118726, 12.75673347, 23.216127…
We can see that i) our priors are relatively weak as the predictions span several orders of magnitide, and ii) our empirical data falls within this range.
Compute Posterior
Fit Model to EEG Data
We will now fit the model to the data.
## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: rms ~ wg - 1 + (1 | subject)
## Data: d_eeg (Number of observations: 400)
## Samples: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
## total post-warmup samples = 20000
##
## Group-Level Effects:
## ~subject (Number of levels: 25)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.38 0.06 0.28 0.52 1.00 1622 3575
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## wgCM -0.48 0.09 -0.66 -0.30 1.00 720 1453
## wgCMM -0.12 0.09 -0.29 0.06 1.00 708 1354
## wgP2 -0.94 0.09 -1.12 -0.77 1.00 704 1308
## wgP3 -0.82 0.09 -0.99 -0.64 1.00 697 1411
## wgP31M -0.43 0.09 -0.61 -0.25 1.00 723 1402
## wgP3M1 -0.14 0.09 -0.31 0.04 1.00 711 1410
## wgP4 -0.48 0.09 -0.65 -0.30 1.00 693 1360
## wgP4G -0.25 0.09 -0.43 -0.07 1.00 700 1346
## wgP4M 0.16 0.09 -0.02 0.33 1.00 710 1450
## wgP6 -0.48 0.09 -0.65 -0.30 1.00 712 1407
## wgP6M -0.04 0.09 -0.22 0.14 1.00 699 1445
## wgPG -0.93 0.09 -1.10 -0.75 1.00 693 1457
## wgPGG -0.72 0.09 -0.89 -0.54 1.00 699 1333
## wgPM -0.59 0.09 -0.76 -0.41 1.00 703 1366
## wgPMG -0.26 0.09 -0.43 -0.08 1.00 687 1320
## wgPMM -0.07 0.09 -0.24 0.11 1.00 699 1382
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.23 0.01 0.21 0.25 1.00 10720 11643
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
We will now look at the model’s predicts for the average participant (i.e, ignoring the random intercepts).
Fit Model to Psychophysical Data
We will now fit the model to the data.
## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: threshold ~ wg - 1 + (1 | subject)
## Data: d_dispthresh (Number of observations: 186)
## Samples: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
## total post-warmup samples = 20000
##
## Group-Level Effects:
## ~subject (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.39 0.11 0.23 0.65 1.00 3033 5669
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## wgCM -0.22 0.17 -0.54 0.12 1.00 2113 3920
## wgCMM -1.15 0.16 -1.47 -0.82 1.00 2105 4043
## wgP2 -0.29 0.17 -0.61 0.06 1.00 2245 4824
## wgP3 -0.69 0.17 -1.00 -0.35 1.00 2079 3574
## wgP31M -1.18 0.17 -1.49 -0.84 1.00 2251 4091
## wgP3M1 -1.35 0.17 -1.67 -1.01 1.00 1993 3771
## wgP4 -0.89 0.16 -1.21 -0.55 1.00 2248 3948
## wgP4G -1.22 0.17 -1.54 -0.89 1.00 2046 4164
## wgP4M -1.29 0.17 -1.60 -0.95 1.00 2191 3908
## wgP6 -1.20 0.17 -1.52 -0.86 1.00 2063 3922
## wgP6M -1.42 0.17 -1.74 -1.07 1.00 2217 3949
## wgPG 0.36 0.17 0.04 0.71 1.00 2020 3794
## wgPGG -0.31 0.17 -0.62 0.03 1.00 2068 3782
## wgPM -0.79 0.17 -1.10 -0.43 1.00 2488 4261
## wgPMG -0.96 0.16 -1.28 -0.63 1.00 2430 3675
## wgPMM -1.17 0.17 -1.48 -0.83 1.00 1882 3745
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.42 0.02 0.37 0.46 1.00 14581 14557
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
EEG Control Data
We will also fit models to the control data. As we can see from Figure 2.4, the group differences are much smaller.
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess